Nucleation barriers in tetrahedral liquids spanning glassy and crystallizing regimes 
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Crystallization and vitrification of tetrahedral liquids are important both from a fundamental 
and a technological point of view. Here, we study via extensive umbrella sampling Monte Carlo 
computer simulations the nucleation barriers for a simple model for tetrahedral patchy particles in 
the regime where open tetrahedral crystal structures (namely cubic and hexagonal diamond and 
their stacking hybrids) are thermodynamically stable. We show that by changing the angular bond 
width, it is possible to move from a glass-forming model to a readily crystallizing model. From the 
shape of the barrier we infer the role of surface tension in the formation of the crystaUine clusters. 
Studying the trends of the nucleation barriers with the temperature and the patch width, we are 
able to identify an optimal value of the patch size that leads to easy nucleation. Finally, we find 
that the nucleation barrier is the same, within our numerical precision, for both diamond crystals 
and for their stacking forms. 



o 



o 
o 



> 

00 
00 

cn 
o 



X 



I. INTRODUCTION 

Crystallization is central to several fields, from materi- 
als research to biological science, not to mention its tech- 
nological relevance [l|, [l] • Several human pathologies are 
also caused by crystal nucleation in protein solutions Q . 
Understanding crystal nucleation requires both the eval- 
uation of the stability fields of the fluid and crystal phases 
(i.e. knowledge of the chemical potentials of all possible 
phases), as well as the evaluation of the thermodynamic 
barriers controlling the formation of the stable phase and 
of the kinetic pre-factors fixing the timescale of the dif- 
fusional processes. 

In recent years, several numerical methodologies have 
been developed for accurately evaluating phase diagrams 
from the free energies of the fluid and crystal phases. We 
refer the reader to the review by Vega and coworkers [1] . 
Also for crystallization, various methods are now avail- 
able for calculating free energy barriers and nucleation 
rates making it possible to generate accurate data 

for model potentials and, more importantly, to compare 
the numerical results with theoretical predictions, mostly 
based on classical nucleation theory (CNT) |10l - [l5| . as 
well as, when possible, with experimental data. 

One of the main motivations for the development of 
the special methods for studying nucleation is the proper 
sampling of the equilibrium cluster size distribution N{n) 
within the metastable liquid, i.e. the number of crystal- 
like clusters composed of n particles. The work of forming 
a cluster of size n from the liquid is given in terms of N[n) 
by i, 



l3AG{n) 



In 



N{n) 



(1) 



where Np is the number of particles in the system and 
/3 = {kBT)~^, where T is the temperature and fee is 
the Boltzmann constant. Below the melting tempera- 
ture, /3AG(n) has a maximum value /SAG* at a critical 



nucleus size n* . (3 AG* enters into the expression for the 
rate of nucleation exponentially, and is therefore the main 
consideration in determining the rate from a thermody- 
namic perspective. Very generally, the larger the differ- 
ence in chemical potential between liquid and crystal A/i, 
often referred to as the driving force for crystallization, 
the smaller (3AG* is. Conversely, the surface tension 
7 between crystallite and surrounding liquid always acts 
against the growth of crystallites, and so (3 AG* increases 
with 7. Hence, the ability of a system to crystallize is 
governed by the interplay between A/x and 7. 

Recently, simulation studies have begun to address 
crystallization of tetrahedral liquids [l6l - l2ll | . This class 
of interesting liquids includes biological and technologi- 
cally important molecular and atomic materials such as 
water, silica, silicon and carbon. One common feature 
in these systems is the formation of open, low density 
structures such as the diamond cubic (DC) crystal. For 
carbon, simulations have shown the importance of the liq- 
uid structure in governing nucleation barriers (22l . [23| . In 
a generalized form of the Stillinger- Weber model for sili- 
con [i^l , the degree of tetrahedrality was shown to have 
a strong impact on DC nucleation and its interrelation 
with the appearance of a metastable liquid-liquid criti- 
cal point |25| . Further, for silicon and germanium, the 
lower density of the DC crystal with respect to the liquid 
was shown to give rise to a preference for nucleation near 
the surface, with similar implications for water [2^. The 
self-assembly of a DC colloidal crystals is also of interest 
in the field of photonics, since such an ordered structure 
of dielectric spheres is expected to exhibit a band struc- 
ture with a large gap within the visible wavelengths of 
light m. 

Tetrahedral liquids are interesting from another per- 
spective as well. Often, in these systems crystal forma- 
tion is inhibited by the onset of a glassy behavior which 
dramatically slows down the microscopic dynamics, thus 
preventing the transition to the ordered structure. The 
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glass forming ability of silica and the ease with which 
water crystallizes stand in stark contrast to each other. 

Motivated by these issues of scientific and technolog- 
ical interest, recently, some of us [l^ have thoroughly 
investigated the question of how to best design tetrahe- 
dral patchy colloidal particles that spontaneously assem- 
ble into open crystal structures through computer sim- 
ulations of the Kern-Frenkel (KF) model [20|. It was 
shown that the angular width of the patch plays a key 
role in determining if the system, at low temperatures, 
forms a disordered glass structure or a crystal, a result 
which may also be of interest for interpreting the glass- 
forming ability of atomic and molecular systems. When 
the patch width is small, the crystal coexists with a fluid 
and /3|A/i| (the driving force for crystallization) increases 
quickly with undercooling, while when the patch width 
is large, the crystal coexists with a fully bonded liquid 
state and /3|A/i| grows only moderately with undercool- 
ing. It was also shown that the model spontaneously 
forms a crystal composed of a mixture of two open tetra- 
hedral structures, the DC and the diamond hexagonal 
(DH) crystals, and that the chemical potentials for the 
two polymorphs were indistinguishable within the preci- 
sion of the calculations. 

In this article we proceed one step further by investi- 
gating, for the same tetrahedral model, the patch width 
dependence of the free energy barrier to nucleating the 
DC and/or DH crystals. Interestingly, we find that, at 
comparable undercooling, the barrier height 13 AG* docs 
not monotonically decrease with increasing patch width 
as one would expect from a consideration of /3| A/z| alone. 
Indeed, comparing the barrier shape with CNT predic- 
tions, we find that the surface tension 7 increases on de- 
creasing patch width, a result that we tentatively connect 
to the larger structural and density difference between 
fiuid and crystal. We do find that for narrow patches 
I3AG* generally shows a strong T dependence, rapidly 
decreasing to a homogenous nuclcation limit, while for 
wide patches f3AG* decreases more slowly with decreas- 
ing T, allowing the system to reach the glassy regime. 
Additionally, at the widest patch studied, the barriers 
remain very large (of the order of 50 fc^T) even for signifi- 
cant undercooling. This confirms that despite the benefit 
of a lower surface tension, the lack of a buildup of a large 
difference in chemical potential is mainly responsible for 
disfavoring crystallization of particles with wide patches. 
The presented results confirm that DC/DH will spon- 
taneously form when the angular width of the patch is 
sufficiently small. 

The remainder of this article is organized as follows. In 
Section II, we outline our criteria for defining crystal-like 
clusters and the umbrella sampling Monte Carlo proce- 
dure for calculating nucleation barriers. In Section HI 
we present our results, including a fairly detailed exami- 
nation of the robustness of (3 AG* to varying the cluster- 
defining criteria and the dependence of n* and cluster 
composition (e.g. percentage of DC particles) on them. 
We then present our summary and conclusions in Section 
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FIG. 1: Equilibrium phase diagrams for all the models stud- 
ied in the low P, low T region, where the open crystal struc- 
tures are stable. The short blue lines correspond to the tem- 
peratures for which the nucleation barriers have been investi- 
gated. The circles mark the gas-liquid critical point. For the 
narrow patch models [(a) and (b)], nucleation is so effective 
that the location of the critical point can not be properly es- 
timated. However, we have checked that the studied isobar is 
located significantly above the critical pressure. 
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II. METHODS 

We perform biased Monte Carlo (MC) simulations [s^l 
at constant T and pressure P of the Kern-Frenkel 
model [2^ with Np — 1000 particles, each with hard 
sphere diameter cr, and each having four tctrahcdrally 
arranged attractive patches of range 5 = 0.24(7, strength 
Mo and angular width 26* [H [H. We report T and P 
as dimensionless quantities, after rescaling by uo/kB and 
uo/a^, respectively. We investigate four different values 
of 9, namely those corresponding to cos 6* = 0.98, 0.96, 
0.94 and 0.92 {6 = 11.5°, 16.3°, 19.9° and 23.1°). For the 
value of S we use, the condition cos 9 > 0.9151 guarantees 
that a patch can only accommodate a single bond, and 
therefore the maximum number of bonds per particle is 
restricted to four 
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Fig. [T] shows the phase diagrams for the Kern-Frenkel 
models we study here, obtained either directly from 
Ref. [2^ [panels (a) and (d)] or by extending the cal- 
culations to cases not reported in that work. The phase 
diagrams show high and low density fully bonded crys- 
tal phases, in which each particle is bonded with four 
neighbors (i.e. the energy per particle is —2uq). The 
low density crystal is an open tetrahedral structure and 
it can exist in two polymorphs, DC and DH. As noted 
in Ref. [2^, the chemical potentials for the DH and DC 
structures are largely indistinguishable, and indeed it has 
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been observed that when the Hquid crystalhzes sponta- 
neously to the open tetrahedral structure, it does so by 
forming a stacking of DC and DH planes [23 • Indeed, DH 
and DC crystalline layers can stack in an analogous way 
to the FCC/HCP stacking in hard-spheres. The dense 
phase, composed of two interpenetrating fully bonded 
DC structures, is a body centered cubic (BCC) crystal. 
The fluid separates into gas and liquid phases below the 
critical temperature. Since the range of the potential 
is short compared to the particle size, the critical point 
is almost always metastable with respect to the crystal 
phase. 

In the following, we study the crystallization barriers 
to the DC/DH crystal at one selected pressure. Specifi- 
cally, we choose P = 0.03 in order to avoid interference 
with gas-liquid phase separation. For all the models stud- 
ied, with the exception of cos 6* = 0.98, the most stable 
phase at that pressure is the open tetrahedral crystal. 
In the case of cos9 = 0.98, BCC is the stable phase at 
this pressure, but as previous studies have shown, spon- 
taneous crystallization results always in the open struc- 
ture dl], an example of the Ostwald step rule [30l - [3^ . 

The nucleation to tetrahedral crwtals has been stud- 
ied previously in simulations [12, [H, [1^ and we follow 
the established methodology to evaluate the free energy 
barriers. A novel issue arises from the presence of two 
crystals with different symmetries in the same stability 
field, as we discuss in the following. 

We use Stcinhardt bond order parameters [l^ based 
on spherical harmonics of order I = S. For each particle 
we define the complex vector 

mm{i) = j^-rrr Yimihj), (2) 

where the sum is over the Ni,{i) neighbors of particle i, 
defined as those particles within a distance of (1 + S)a = 
1.24(7 of particle i. The dot product 

Cy = ^;m(i)9L(j)) (3) 

m.— — l 

where 

/ ' 

qhn{i) = qhn{i)/ { Y l9'™(*)n (4) 

\m=-l ) 

and (?;*„(i) is its complex conjugate, determines the de- 
gree of orientational correlation between neighboring par- 
ticles i and j. Fig. [2] shows the distributions for the 
fluid and the DC and DH crystals. The DC distribution 
is peaked around = — I only, while the DH crystal 
also shows a peak near Cy = —0.1. In the DH crystal, 
each particle has three neighbors with Cy « — 1 and one 
with Cij ~ —0.1, while in the DC crystal all four bonded 
neighbors have cij w — 1. This provides a local basis for 
distinguishing particles as being DC or DH. The fluid dis- 
tributions are very wide and show sharp peaks for large 



values of cos B (small bonding angles) . The peaks become 
more intense on heating, which we attribute to a lower 
density and higher energy; we identify the peaks as sig- 
nals of specific geometrical assemblies in the fluid with 
unfilled bonds. For example, a dimer has = — 1, while 
both bonds in a trimer have — —0.82 

The possibility of separating the ranges in which 
crystal-like or fluid-like particles are mostly contributing 
offers a way of associating a value of Cy- with a local 
structure (crystal or fluid). Usually a threshold number 
of crystal-like connections is selected to distinguish par- 
ticles as being fluid-like or crystal-like. In the following, 
we start by dcflning a crystal-like connection between 
neighbors as one with < q^, with — —0.87 and 
a crystal-like particle as one which has three or more 
crystal-like connections. This definition does not allow 
one to discriminate between DC and DH and hence ap- 
pears to be a reasonable first step in the investigation of 
the nucleation barriers [22, 23, 2^. We complement this 
study with an additional investigation where we differen- 
tiate between DC and DH by requiring that a solid-like 
particle have four neighbors with Cy < —0.87 in the DC 
case or three neighbors with Cy < —0.87 and one with 
—0.3 < Cij < 0.1 in the DH case. For completeness, we 
probe three cases: (i) the case where the growing crystal 
is composed only of DC particles; (ii) the case where the 
growing crystal is composed only of DH particles; (iii) 
the case where the growing crystal is composed of DC or 
DH particles. 

We follow the standard methodology for dcflning a bi- 
asing potential which helps the formation of crystalline 
clusters. To this aim we add to the Kern-Frenkel poten- 
tial a perturbation given by 

= K(nmax - J^o)^, (5) 

where k is a suitably chosen constant that controls the 
range of sampled crystal cluster sizes, centered near ng, 
and rtmax is the size of the largest crystalline cluster in 
the system. Depending on the steepness of AG{n), we 
adjust the value of k to obtain good sampling, but for 
almost all simulations, k = 0.075. To dcflne a crystal- 
like cluster, we state that two neighboring crystal-like 
particles are part of the same cluster. 

New configurations in the umbrella sampling (US) MC 
chains are generated as follows. Given a starting configu- 
ration with largest cluster n^axi perform a trajectory 
of 20 Metropolis MC steps, where one such step is on av- 
erage A^p — 1 particle translations and rotations and one 
volume change, in order to arrive at a configuration with 
largest cluster nl^^^. The new configuration is accepted 
with probability min(l, exp {-^[^(nJn^^) - 0(71^^^)]}). 
If the new conflguration is accepted, it becomes the start- 
ing point for the next MC trajectory. If it is rejected, the 
entire trajectory is discarded and the n^j^^ configuration 
is kept and used as the starting point for generating an- 
other MC trajectory. For the slowest state points, we 
perform 10^ US attempts. We perform several US sim- 
ulations for different values of ng. Generally, two adja- 
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FIG. 2: Probability distributions of the dot product Cij in the 
liquid and crystal structures. The distribution for DC, shown 
only in panel (a), is a single peak near -1 and is nearly indis- 
tinguishable from the DH peak at -1 on the scale of the plots. 
The distributions for the crystals do not vary significantly 
over our T range. 



cent US windows are spaced by Ang = 3 to ensure good 
sampling of the reaction path. For each US window we 
then evaluate the solid cluster size distribution N{n) and 
Pmax('^), the probability that the largest cluster in the 
system is of size n. 

The cluster size distribution N{n) in the NPT ensem- 
ble for each window is worked back from the biased en- 
semble as in Ref. M, with 



N{n) = (cxp [/?(^(n,„ax)]^(n) 



biased 



(6) 



up to a multiplicative constant, where (...) indicate an 
ensemble average. Portions of f3AG{n) are determined 
from simulations at different values of uq up to additive 
constants, in agreement with Eq. [T] These pieces are 
matched by minimizing the difference between overlap- 
ping portions after discarding data for which Pmax('T-) is 
less than 0.01 (to ensure good sampling). Alternatively, 
we find the free energy difference AG{n) — AG{n — l) as a 
weighted average of the values obtained with simulations 
at different uq in which clusters of size n and n ~ 1 have 
been sampled, where the weight is given by cxp[— 2fc(ri — 
no)]. As a result AG(n) = ELi[AG(j) - AG(j - 1)]. 
The two procedures give equivalent results. We then shift 
the curves so that AG(0) = 0. While this is not formally 
correct, since the number of liquid-like particles [A^(0)] is 
not precisely equal to A^p, the error is negligible (of the 
order of 0.01 /cbT or less). We stress here that at this 
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FIG. 3: Nucleation barriers /3AG(n) for various bond angular 
widths and temperatures. Panel a) cos 9 = 0.98, from bottom 
to top T = 0.140, 0.142, 0.143, 0.144, 0.145, 0.146, 0.148, 
0.150. Panel b) cosO = 0.96, from bottom to top T = 0.154, 
0.156, 0.158, 0.160, 0.161, 0.162, 0.163, 0.164, 0.165. Panel 
c) cos6 = 0.94, from bottom to top T = 0.154, 0.156, 0.158, 
0.160, 0.162, 0.164, 0.166, 0.167, 0.168. Panel d) cos9 = 0.92, 
from bottom to top T = 0.154, 0.156, 0.158, 0.160, 0.162, 
0.164, 0.166, 0.168. 



stage we do not differentiate between DH and DC parti- 
cles. We have checked that indeed the largest crystalline 
cluster is a mixture of the two crystal local structures. 
We will address this point in more detail later on. 

To ensure independence of the initial conditions, we 
construct our barriers starting from two different sets 
of data. In the first set we use the biasing potential 
to grow clusters in the fiuid, starting the US simulation 
at no from a configuration extracted from the previous 
windows. In the second set, we seed the biased simu- 
lations at cos6' = 0.98 with crystallite-containing con- 
figurations from homogeneously nucleating unbiased MC 
runs at cos9 = 0.98. Simulations at smaller values of 
cos9 are started from the cos9 — 0.98 data set, after 
equilibration at higher values of T. 



III. RESULTS 

A. Barrier profiles 

In Fig.[3]we show our results for the AG(n) profiles for 
the four models of different patch widths at various T. 
The range of T where the barrier can be calculated with 
the present methodology is restricted both from above 
and below for different reasons. At large T (small su- 
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percooling) the critical nucleus is large compared to the 
system size. At low T, different reasons conspire against 
a proper evaluation of the barrier height. For example, 
in the case cos^ = 0.92, the slow kinetics associated with 
the proximity of a glass transition prevent proper equi- 
libration already at temperatures where the critical nu- 
cleus is of comparable size to the studied system. In this 
same case, it has been suggested on the basis of the ex- 
plicit evaluation of the difference in chemical potential 
between the crystal and the fluid that, at odds with the 
standard behavior, barriers do not grow with further su- 
percooling. The difficulty of evaluating barriers for the 
wide patch model is consistent with this proposition. In 
the case of narrow patches, cosO = 0.98, we are limited 
to a barrier height of the order of 30 /cbT- For lower T, 
despite slow dynamics not being an issue, the results of 
the calculations become more and more affected by the 
presence of secondary crystal clusters. While in theory 
the presence of more than one cluster is not a problem, in 
our opinion this effect highlights an incorrect definition 
of crystallinity at the local level, which artificially breaks 
a single cluster into two (or more) pieces. Hence, at the 
lowest T, the barrier calculations become difficult owing 
to the appearance of clusters of comparable size to the 
largest one in the system. A visual inspection of con- 
figurations confirms the formation of secondary clusters 
occurring next to the primary one, separated by particles 
that fall outside the cutoffs defining crystal-like particles. 
We recall that for a proper evaluation of the barriers in 
the low T limit, where homogeneous nuclcation is tak- 
ing place, one can apply the methodology based on mean 
first-passage times 3^- iO] ■ Finally, we note that we have 
checked for state points where the barrier is of the order 
of 20 k^T the quality of our calculations by comparing 
the US results with N{n) evaluated in unconstrained sim- 
ulations. 



B. Fits 

Within the phenomenological framework of CNT, the 
work of forming an n-sized cluster can be written as. 



l3AG{n) ^ -l3\An\n + /S-fA, 



(7) 



where A is the surface area of the cluster. We therefore 
fit our profiles to 



(8) 



where a = /3\Afi\ and b ^ (3^ and we have assumed a 
compact cluster. Below, we compare our results for a 
against /3|A/xJ, as some studies suggest a very good cor- 
respondence 0,11,1411, and justify our assumption that b 
is equal to (3^ up to a constant factor which depends on 
the shape and density of the crystallites. Recently, it has 
been shown for a soft-core colloidal model that the pro- 
cedure we follow here for determining n may not be suf- 
ficient for describing crystal-like structures, but that the 




FIG. 4: Fit of two of the barrier profiles studied to the CNT 
form (Eq.[8]). While tliere is a deviation at small n, the overall 
representation of the curve is fairly good, especially in terms 
of the barrier height AG* and size of the critical nucleus n* . 
Inset: low-n region of the same data, hilighting the deviation 
from the CNT form. 



scaling represented by Eq.[S]may still hold, i.e., changing 
the parameters used to define a crystalline cluster will 
change the numerical values of a and b [i^l . We explore 
this in some detail below. 

We show two representative /3AG{n) curves along fits 
to Eq.[H]in Fig.U] While there is a deviation between data 
and fit at small n, the overall description is rather good. 
Moreover, the values for n* and AG* extracted from the 
fits coincide with those obtained directly from the bar- 
riers. Improvement of the small n description could be 
accomplished along the lines suggested in Ref. [4^ but it 
is outside the scope of the present work. 

The fits allow us to plausibly extrapolate our AG{n) 
profiles to values of n larger than what we can simulate 
in our A^p = 1000 system, allowing us to estimate n* and 
(3 AG* for all the T that we study. While we do not rely 
strongly on the accuracy of the extrapolations, they do 
provide a stronger sense of the trends in the data. 



C. T dependence of AG* and n* 

Fig. m^a) shows /3AG* as a function of T/T„i, where 
Tm is the melting temperature for the four models [i^ 
[(008 6*, r„i): (0.98, 0.153), (0.96, 0.169), (0.94, 0.172), 
(0.92, 0.174)]. For small angles, the T dependence of 
/3AG* is significant, and can be modeled rather well in 
the present range of barriers with an exponential func- 
tion. The fast decrease of /3AG* provides evidence for 
the inevitability of crystallization. For larger angles, we 
observe both a significant change in slope, i.e., a much 
slower decrease of the barrier height with supercooling, 
as well as a progressive increase of the barrier height at 
fixed supercooling with increasing angle. The trend is ac- 
counted for by the results reported in namely that 
/3| A^l becomes a weaker function of temperature as 6 in- 
creases. On the basis of the results presented in Fig.O^a), 
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FIG. 5: (a) Barrier height /3AG* as a function of T/Tm. For 
high values of cos^, the barrier decreases significantly with 
decreasing temperature, while the slope is much more shal- 
low for low values of cos 9. The two points at lowest T for 
cos© = 0.98 are obtained from US simulations employing a 
more relaxed dj cutoff of q„ = —0.65 after checking for consis- 
tency at higher T. (b) Critical cluster size n* as a function of 
the reduced temperature. Also in this case the trend with T is 
stronger in narrow-patch models than in wide-patch models. 



the model with cos6' = 0.96 appears to be the optimal 
candidate for crystallization from a thermodynamic per- 
spective, since modest supercooling is sufficient to induce 
barriers of height of the order of 10 fc^T. 

Fig. Eljb) shows n* as a function of T/Tm, showing 
similar trends as for I3AG*{T). The cos 6* = 0.98 and 
0.96 models are similar, attaining small n* for a relatively 
small degree of supercooling, while the models with wider 
patches seem to require larger sizes of critical nuclei as 
T decreases. The T dependence of n* is much flatter, 
suggesting that the critical nucleus remains significantly 
large even under deep supercooling. 
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FIG. 6: (a) Barrier height as a function of n^^'"^ (proportional 
to the surface area of the cluster). All models show a linear 
behaviour, (b) Surface coefficient b of the CNT fit as a func- 
tion of the reduced temperature. The surface term does not 
change significantly with T, as compared to the barrier height 
[see Fig. [5ja)] . (c) Reduced number density pa^ as a function 
of T at P = 0.03 for the models studied. Large open symbols 
indicate the range of state points used in the calculation of b. 
Curves are a guide to the eye. 



D. Surface tension 

To illustrate the dependence of the surface tension on 
patch width, motivated by the relation 3/3 AG* = bn*^^^ 
that follows from Eq. |51 we plot in Fig. f3AG* as 
a function of n*'^/^. For all models, a reasonable linear 
dependence is observed, with a slope that progressively 



increases on with cos 6* beyond 0.94. This results in a 
larger critical size for a given barrier height as the patch 
width increases. The near linearity also suggests that the 
surface tension does not have a strong dependence on su- 
percooling. In Fig.ini^b) we show 6 as a function of T/T„i 
as obtained from the fitting procedure for the different 
state points studied. The surface tension increases sig- 
nificantly on decreasing the angular patch width. Once 
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FIG. 7: CNT bulk fit parameter a as a function of temper- 
ature (points) compared witli the difference in ciiemical po- 
tential /3 |A/^| (lines). While the two quantities show similar 
trends, to achieve quantitative agreement we multiply values 
of a in the plot by a universal scaling factor of 0.4. 

again, however, the two models with the largest patches 
are quite similar in their behavior. 

It is worthwhile noting that the densities at our simu- 
lated P for these two wide patch models (at our lowest 
T) are very similar with values near pa^ = 0.47, almost 
matching the crystal density (where p is the number den- 
sity). The density is smaller for cos = 0.96, and smallest 
for cosO = 0.98. Thus, the values of b for the different 
models appear to reflect the differences in density, with 
a better match in density between fluid and crystal giv- 
ing rise to a smaller surface tension. However, a plot 
of the T dependence of the densities in Fig. |6l^c) shows 
that density does not solely determine h, as indicated by 
the variation of density without a correspondingly large 
variation in b within each model. 



E. a vs /3A/i 

In Fig. [7] we plot /?|A/x| [i^ and a (obtained from the 
CNT fits) as functions of T . For comparison purposes, 
all the values of a appearing in the figure have been mul- 
tiplied by a factor of 0.4, i.e. the values of a are signif- 
icantly higher than expected from independent calcula- 
tion of /3|A/i|. However, it is quite comforting that the 
rescaled a matches /3A/i quite well across T and 9, as 
this strengthens our assumption that the fit parameter b 
is proportional to /37 with the same proportionality con- 
stant across all the models. 

On the other hand, it is somewhat perplexing that such 
a large rescaling is required at all. In a few previous stud- 
ies, comparison between a and /3A/i have shown that in 
some cases a close correspondence is observed @, |4l| , 
while in principle the value of a should vary with the 
definition crystal-like particles. There are several po- 
tential reasons why such a disagreement can be found. 
First of all, the assumption of CNT may not be fully 
satisfied, including those concerning the structure of the 
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FIG. 8: Comparison of the nucleation barriers obtained with 
two different order parameters, = —0.65 (dashed lines) 
and g„ = —0.87 for the cos^ = 0.98 model. As expected, the 
barrier heights are consistent (see Ref. [iJl) but the critical 
cluster sizes differ significantly. The more relaxed order pa- 
rameter leads to a larger critical cluster, possibly due to the 
fact that it includes more particles on the surface. In the inset 
we also show the dependence of the bulk CNT term a as a 
function of the order parameter, showing that it increases by 
making the order parameters more strict. 

crystallites and the sharpness of the interface [l^] ■ Sec- 
ond, an imperfect classification of particles as solid-like 
leads to an imperfect evaluation of the cluster size and a 
classification-dependence of a and h. 

As a test of this second hypothesis and to illustrate 
the effect of the parameters chosen to define crystalline 
clusters, we plot in Fig. [S] sets of barrier profiles obtained 
for cos Q ~ 0.98 using g„ = —0.65 compared against using 

= —0.87. While AG* remains invariant to within er- 
ror with respect to at all investigated T, n* increases, 
as expected, for smaller values of qu- In other words, the 
size of the critical nucleus is strongly dependent on the 
criteria chosen to define solid-like particles. It is not a 
surprise that AG* is a much more robust value, since the 
work required to produce the critical nucleus is indepen- 
dent on how the cluster is described. In other words, the 
"real" critical nucleus (i.e. configurations which 50 per 
cent of the time will crystallize and 50 per cent will melt 
again into a fluid state) is what needs to be sampled in the 
simulation. The number of particles that compose this 
cluster depends on the definition, affecting the perceived 
size, but not the barrier height. These conclusions are in 
agreement with the work of Filion et al [4^ on spherical 
particles. 

In the inset of Fig. [51 we show a as a function of 
for T = 0.146 and cos 9 = 0.98 and see that a varies by 
roughly a factor of 2 for the reasonable range of g„ we 
have explored. Coincidentally, the value for /3|A/i| for 
this state point is 0.56. Thus, for an appropriate choice 
of Qu, Eqs. [7] and [5] describe the nucleation barrier profiles 
both in scaling and in the numerical value of Afi, but how 
to choose the appropriate g„ a priori is not clear. 

The more inclusive choice of g„ = —0.65 allows one to 
more easily explore the regime of low barrier heights, less 
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FIG. 9: Schematic description of the criteria that require 
four connections for defining solid-like particles. The values 
of Cij are displayed next to the bonds that the two central 
particles form. The particle on the left has three bonds with 
Cij < —0.87 and one bond with —0.3 < Cij < 0.1, registering it 
as a DH particle. The particle on the right has four bonds with 
Cij < —0.87, registering it as a DC particle. In the mixed four- 
connections case, both particles would be considered solid- 
like. 
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FIG. 10: Nucleation barriers for pure DC, pure DH, and 
mixtures of the two. For T — 0.144 the four-connections 
curves (pure DC, pure DH and mixed DH-DC) are all quite 
similar to each other and yield the same barrier height as 
the three-connections case, but smaller n* . Barriers at T = 
0.142 compare the mixed three- and four-connections cases 
and confirm the invariance of AG* and the change in n* with 
different criteria for defining clusters. 

hampered by the formation of secondary clusters. Using 
Qu = —0.65 we are able to add two additional points to 
the curve for cos6 ~ 0.98 in Fig. [SJa). 



F. DH vs DC 

As we have alluded to before, the clusters are formed 
by a mixture of local DC and DH particles. Indeed, 
the criterion of having at least three connected neigh- 
bors {cij < Qu) for defining solid-like particles, accepts 



either a locally DC or DH particle as solid. However, this 
definition is more likely to misidentify a DH-like particle 
as liquid-like, since particles in the DC crystal have four 
such connections, while DH particles only have three. 
The biasing potential may therefore tend to grow clusters 
richer in DC than would occur naturally. In the follow- 
ing we therefore compare this criterion with other criteria 
which enforce the growth of pure DC, pure DH or an un- 
biased mixture of the two. We do this at cos6' = 0.98, 
where the calculations are the fastest. 

For pure DC, we retain = —0.87 and simply re- 
quire that the number of connections be exactly four (as 
opposed to at least three). With this more restrictive 
criterion, the solid-like particles must have a DC envi- 
ronment. For pure DH, we require exactly three connec- 
tions with Cij < qu = —0.87 and one connection with 
—0.3 < Cij < 0.1. In this way, DC particles are excluded 
and can not participate in crystalline clusters. For the 
mixed DH/DC case, a particle is solid-like if it is deter- 
mined again by exactly four connections, but this time 
these four connection must satisfy either the criterion for 
a DC particle or the criterion for a DH particle. The 
cartoon in Fig. [HI provides a graphical explanation of the 
criteria for identifying solid-like particles. 

In Fig. [ini we show the barrier profile results at 
T = 0.144 with cos 6* = 0.98 for the four different defini- 
tions of solid-like particles (three-connections, pure DC, 
pure DH and unbiased four-connections mixture), in or- 
der to asses if the energy barriers of the pure crystals 
are different from that of the mixture and to ascertain 
the difference between using three or four connections in 
the mixed case. We can not discern any real difference 
in terms of the barrier height between the four criteria. 
This interesting observation supports the view for this 
short-range model, that in addition to the DH and DC 
free energies being essentially identical, the surface ten- 
sions are also similar. Even more, if we restrict ourselves 
to the comparison of all criteria in which four connec- 
tions are required, then not only is the barrier height 
identical, but also (within our numerical precision) the n 
dependence of the barrier profile. This suggests that the 
pathways to crystallizing DC and DH are quite similar 
as well. The pure DH case appears to produce a slightly 
larger n* , but this will also likely depend on the bounds 
on Cij near —0.1. Comparing the three-connections cri- 
terion with the more restrictive four-connections cases, 
we sec that the four-connections criterion simply results 
in an apparently smaller critical cluster size, providing 
one additional piece of evidence for the sensitivity of the 
profile on the solid-like definition accompanied by the in- 
sensitivity of the barrier height. 

We repeat the comparison at the lower T = 0.142. 
However, we are unable to construct barriers for the pure 
DC and DH cases. At this temperature and even at mod- 
erate values of uq, the presence of a pure DC cluster ap- 
pears to act as a template for DH growth (and viceversa) . 
Essentially, the DH particles, which are not counted as 
solid-like, become part of the growing crystal, providing 
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a bridge between only apparently different DC crystal 
clusters. Thus, the system crystallizes while still regis- 
tering a small largest cluster in the system and our order 
parameter no longer describes nucleation. 

This issue is greatly reduced for the four-connection 
mixed case and we are able to calculate the barrier at 
T = 0.142. Again, the barrier height is comparable to 
the one calculated for the case of three connections. At 
lower T, also for the four-connections mixed case, it be- 
comes impossible to properly evaluate the barrier, since 
the overly restrictive criterion will misidentify as liquid- 
like particles sufficiently crystalline to participate in crys- 
tal growth. 

While the barrier heights are the same across the dif- 
ferent cluster criteria, the DC/DH composition of the 
clusters varies. To quantify this in a basic way, we ana- 
lyze an ensemble of largest clusters extracted from a set 
of US simulations employing the mixed four-connections 
solid-like criterion in a wide range of values. We 
find that 54% of the solid particles are DC (standard 
deviation 17), which is the same as the composition for 
clusters appearin g in spontaneous nucleation studies as 
reported in Ref. [28|. Performing the same evaluation 
starting from US simulations, employing this time the 
mixed three-connections solid- like criterion (g„ = —0.87), 
the fraction of DC particles within the largest cluster is 
73% (standard deviation 22). The result is the same for 
qu = —0.65. This confirms that the choice of at least 
three connections with Cij close to minus one does in- 
deed introduce a bias toward the DC structure. For this 
model, however, this bias does not measurably affect the 
nucleation barrier. 



G. BCC 

For all the patch widths considered except cos = 0.98 
we are working in the stability field for DC/DH. For the 
case of cos 6' = 0.98, the stable crystal phase at P = 0.03 
is the BCC (see Fig. [T]), even if spontaneous nucleation 
indicates that the crystal that forms is the DC/DH mix- 
ture. To verify that the barrier for nucleating the BCC 
is significantly larger than the one for nucleating DC or 
DH or their mixture for cos = 0.98, we evaluate the bar- 
rier at T = 0.142 using the same procedure described for 
the tetrahedral crystals. The definition of solid-like BCC 
particle is based on the / = 6 spherical harmonics Q, 
defining neighbors to be connected when the scalar prod- 
uct of the / = 6 harmonics is greater than 0.5. A particle 
is classified as a solid-like if it has at least six connections. 
We note that Ref. Q does not use normalized qim{i) vec- 
tors for calculating Cij, so we have determined cutoffs 
based on our distributions of Cy and of the subsequent 
number of connections a particle has in the liquid and 
in the crystal. We do find that the barrier to nucleating 
BCC is at least 70 fceT larger than for tetrahedral crys- 
tals, supporting the lack of observation of spontaneous 
BCC nucleation at the pressure we consider here. 
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FIG. 11: Estimates of homogeneous nucleation temperature 
Tx" (crosses) and (circles) as functions of cosO. The 
solid vertical line indicates the value of cos^ = 0.9151 at 
which more than one bond per patch becomes possible for 
S — 0.24. The upper horizontal dashed line marks the ap- 
proximate value of Tg/Tin = 0.8 for silica, while the lower one 
indicates the general value of Tg = 2Tm/?>- Inset shows the 
determination of Ti^ and from the crossing of exponen- 
tial fits to the four lowest T points for each value of cos 9 with 
l3AG* = 10 and 20, respectively. 



IV. SUMMARY AND CONCLUSIONS 

Previous work for this tetrahedral patchy model 
showed that the driving force for nucleation (quantified 
as difference in the crystal and fluid chemical potentials) 
decreases as the patch width increases. Fig. [7] summa- 
rizes this result by showing /3|A^| increasing rapidly to 
large values as T decreases below for narrow patches, 
while increasing only slowly and remaining fairly small 
for wide patches. 

We find that despite this reduced driving force, the 
barriers to nucleation actually are smallest for the cos 9 = 
0.96 case [as shown in Fig. [S^a)], in the sense that 
this model reaches the homogeneous nucleation limit, 
marked for example by the barrier reaching a value of 
/3AG* = 10, at the highest T/T„ amongst the stud- 
ied models. The increased similarity between liquid and 
crystal in terms of energy and density as patch width in- 
creases not only brings the chemical potentials of liquid 
and crystal closer in value (tending to increase the nucle- 
ation barrier), but also reduces the surface tension (tend- 
ing to lower the barrier). Thus, in the range of narrow 
angles where crystallization is readily observed, compe- 
tition between |A/i| and 7 leads to optimal nucleation at 
an intermediate patch width. 

Increasing the patch width beyond cos 6* = 0.94 no 
longer significantly reduces 7, while | A/i| continues to de- 
crease, causing an increase in nucleation barrier heights. 
For cos6' = 0.92, this resulting increase is quite large, 
with (3 AG* estimated to be of order 50 at the lowest T 
that we can simulate. Moreover, the rate of decrease of 
PAG* with T appears to be quite slow, requiring signif- 
icant supercooling to reach accessible nucleation barrier 
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heights. 

The evaluation of the barriers requires a definition of 
sohd-like particles. We have checked that the barrier 
height is essentially insensitive to the exact choice of 
the cutoff used to define solid-like connections, consis- 
tent with the observation in Ref. [31 for hard-spheres. 
The size of the critical nucleus is instead significantly de- 
pendent on the definition of solid-like particles, again in 
agreement with the conclusions in Ref. 01- We have 
also compared several definitions of solid-like particles to 
probe the crystallization of a pure DC, a pure DH and a 
mixed DC-DH structure. Interestingly, we find the same 
value of (3 AG* for all structures and, in this case, similar 
n* . The fact that the entropy gain in creating a mixed 
structure does not appreciable lower the work of forming 
a critical nucleus, but is sufficient to generate a preva- 
lence of mixed structures when spontaneous nuclcation 
takes place [28j . warrants further investigation. 

At this point, we would be remiss if we were not to 
comment on the interplay between dynamics and ther- 
modynamics in controlling the rate of nucleation, and 
indeed governing the glass-forming abilities of the sys- 
tem. We already see an example of this within our data. 
The barriers at the lowest T studied for cos 6* = 0.98 and 
cos 6 = 0.94 are the same within error, and have a value of 
/3AG* ~ 17. However, as mentioned previously, the dy- 
namics (in terms of the diffusion coefficient) are 40 times 
slower for the wider patch case. For simulations, this is 
a significant number. The slow decrease in /3AG* with 
T, combined with the expected slow-down in dynamics, 
suggest that a search for spontaneous nucleation in un- 
constrained simulations would target T not far from the 
lowest T for which we have calculated the barrier. The 
case of cos 9 = 0.92 is more obvious since its dynamics 
are even slower at a given T than those of cos 9 = 0.94: 
there is no hope of seeing nucleation in this model, which 
for all intents and purposes is a glass-former. 

To place these results in slightly broader context, that 
may be illuminating for molecular tetrahedral liquids 



that are quite different in their glass-forming and crys- 
tallizing properties, we use the nearly exponential be- 
havior of (3 AG* to find estimates for the temperature at 
which the homogeneous nucleation limit is reached, Tx- 
We define two such estimates using /3AG{T^) = 10 and 
f3AG{T^) = 20, and determine them by fitting the low- 
est four points in T from Fig. EJa) for each model to 
exponentials and extrapolate where necessary. In Fig. 1111 
we present the resulting estimates of Tx alongside a verti- 
cal line showing the widest patch that guarantees a single 
bond per patch, i.e. four bonds per particle. Also shown 
are two horizontal lines indicating the approximate value 
of Tg/Tm = 0.80 for silica ^ (where Tg is the glass tran- 
sition temperature) , as well as the general rule of thumb 
Tg/T^ = 2/3 [13] which has been shown to be valid for a 
large class of molecular [i^ and polymeric systems [49| . 
For Tx/Tm below Tg/Tm = 2/3, a model may be consid- 
ered a glass-former; our estimates suggest that we may 
be approaching this regime with our widest patch model. 
The presence of a maximum in Tx confirms that decreas- 
ing the angular range for bonding favors glass formation. 
Tetrahedral colloidal particles will thus form crystals only 
if the bonding angular width is small. 

As a final remark, we point out that our findings 
may shed some light on the glass-forming and crystalliz- 
ing abilities of molecular or atomic tetrahedral network- 
forming liquids, contributing insight as to why, for exam- 
ple, water crystallizes while silica more readily forms a 
glass. 
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